home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1994 August / August CD.bin / Shareware / Education / numericalmethods Folder / chap_9 / a9_4.m < prev    next >
Encoding:
Text File  |  1994-06-05  |  1.6 KB  |  69 lines  |  [MATS/MATL]
  1. echo off;
  2. % NUMERICAL METHODS: MATLAB Programs, (c) John H. Mathews 1994
  3. % To accompany the text:
  4. % NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992
  5. % Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A.
  6. % This free software is complements of the author.
  7.  
  8. % Algorithm 9.4 (Runge-Kutta Method of Order 4).
  9. % Section    9.5,    Runge-Kutta Methods, Page 460
  10. echo on; clc; format long; hold off; clear
  11.  
  12. % This program implements the Runge-Kutta method
  13.  
  14. % for solving the initial value problem
  15.  
  16. %     y' = f(t,y)    with    y(a) = y0
  17.  
  18. %    Define and store the function f(t,y)    in the M-file  f.m 
  19. % function z = f(t,y)
  20. % z = (t-y)/2;
  21.  
  22. delete f.m
  23. diary f.m; disp('function z = f(t,y)');...
  24.            disp('z = (t-y)/2;');...
  25. diary off;
  26.  
  27. f(0,0); % Remark. f.m and rk4.m are used for Algorithm 9.4
  28. pause      % Press any key to continue.
  29.  
  30. clc;
  31. %    Place the endpoints of [a,b] in  a  and  b.
  32.  
  33. %    Place the initial value ya = y(a) in  ya.
  34.  
  35. %    Place the number of subintervals in  m.
  36.  
  37. a  = 0;
  38. b  = 3;
  39. ya = 1;
  40. m  = 12;
  41. [T,Y] = rk4('f',a,b,ya,m);
  42. points = [T;Y];
  43.  
  44. pause    % Press any key to see the list of solution points.
  45.  
  46. clc;, clg;
  47. c = 0;
  48. d = 1.75;
  49. axis([a b c d]);...
  50. plot(T,Y,'g');...
  51. hold on;...
  52. plot([a b],[0 0],'b',[0 0],[c d],'b');...
  53. if m<=32,
  54.   plot(T,Y,'or');
  55. end;...
  56. xlabel('t');...
  57. ylabel('y');...
  58. title('Runge-Kutta solution to y` = f(t,y)');...
  59. grid;...
  60. axis;...
  61. hold off;...
  62. shg; pause    % Press any key to continue.
  63.  
  64. Mx1 = 'Runge-Kutta solution to y` = f(t,y).';
  65. Mx2 = '     t(k)               y(k)';
  66. clc,echo off,diary output,...
  67. disp(''),disp(Mx1),...
  68. disp(''),disp(Mx2),disp(points'),diary off,echo on
  69.